Linux shell trick for bioinformatics 系列文章之一

    由一个纯生物背景的人，转到生物信息学领域之后。刚开始学会写脚本的时候，仿佛开启了一扇通往新世界的大门。好多以前觉得很难完成的事，写个脚本轻松完成。随着时间的推移，攒下了越来越多的小脚本。伴随着技能树的拓展，知道了很多事根本不用写脚本，直接用shell一行搞定。

    并且用shell来处理这些小事的时候，还有一个好处就是通过管道符将这些linux操作串联以来，形成一套pipeline。 诚然，写脚本也可以，但是用linux代码更加符合自己对简单，优雅的审美要求，所以告别那些无用的dirty code，拥抱简洁强大的Shell。

1.   将一个文件按行倒序排列，即将第一行变为倒数第一行，第二行变为倒数第二行，一次类推。我们有一个稍微笨一些的方法：

awk 'BEGIN {x=0}{x=x+1}{print x,$0}'  yourfile.txt | sort -nr -k 1| sed 's/^.* //g' |le 

或者我们采用一个优雅至极的方法；

tac  yourfile.txt

就是linux常用命令cat反过来就是倒序输出。

2.  对于Rfam数据库下下来的noncoding RNAs的fasta文件，注释信息就在fasta文件的header里面，如果我们想要提取所有的sequence的注释信息。

grep "^>" rfam.fasta | tr -d ">" |le

３.  用Trinity坐组装的之后得到的fasta格式的转录组文件，我们要做注释的话，有些header可能过长，所以我们要将header给截短些，仅保留最关键的信息。

cut -f 1 -d " " trinity.fasta |le

４.  统计fastq的所有碱基的数目。

awk 'BEGIN{sum=0;}{if(NR%4==2){sum+=length($0);}}END{print sum;}' your.fq

5. SAM文件中根据某一行的数值进行筛选，并且保留header.  以拟南芥这种模式植物为例，和拟南芥的记忆组比对得到的SAM文件有9行header，我们想保留header，再根据某行值进行过滤。

awk '$3<1000 || NR <= 9' your.SAM

6.  将fastq格式文件转化为fasta格式。

sed '/^@/!d;s//>/;N' your.fastq> your.fasta

7.  当我们得到基因的注释之后，得到了two column的基因注释table, 但是里面有很多NA值。但我们不想看到NA，怎么办呢。方法很多。

perl -ne 'print unlesss /NA/'  file.txt

awk  '{if(!/NA/){print $0}}'  file.txt

awk '$0 !~ /NA/ {print}'  file.txt

grep -v "NA"  file.txt

sed -e '/*NA*/d'  file.txt

8.  递归地查找当前目录下所有以"fasta"为后缀的文件，统计其行数。如果想统计fasta文件的所有碱基数目，请参考统计fastq碱基数目那一条。

find . -name "*.fasta " | xargs wc -l

9.  当你有两个gene list的文件，要找出仅在gene.file2中存在的行，不在gene.file1中出现的gene用来做GO分析。

grep -Fxv -f gene.file1 gene.file2

10.  在某个目录下有很多文件，你想看看你最感兴趣的基因名字出现在哪个文件里，文件很多，并且子文件夹还有很多文件，即在一个目录下递归地查找匹配某字符串的文件。怎么办呢。

grep -Hrn "shengxinyuan" .